The design of Kinetic Functionals for Many-Body Electron Systems : Combining 
analytical theory with Monte Carlo sampling of electronic configurations. 
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In a previous work [L. Delle Site, J.Phys.A 40, 2787 (2007)] the derivation of an analytic expression 
for the kinetic functional of a many-body electron system has been proposed. Though analytical, 
the formula is still non local (multidimensional) and thus not ideal for numerical applications. In 
this work, by treating the test case of a uniform gas of interacting spinless electrons, we propose a 
computational protocol which combines the previous analytic results with the Monte Carlo (MC) 
sampling of electronic configurations in space. This, we show, leads to an internally consistent 
scheme to design well founded local kinetic functionals. 
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The design of electronic kinetic density functional is 
becoming a field of rapidly growing interest due to the 
central role played by this latter in Orbital-Free Density 
Functional (OFDFT) methods. Such methods allow for 
large computational advantages compared to standard 
Kohn-Sham (KS) based techniques since they are com- 
putationally much less demanding 0, Q, H, 0, S] ■ More- 
over, as the calculations are performed only in real space 
they allow for an optimal design of an interface with clas- 
sical codes to the aim of connecting the electronic and 
the larger atomistic scale within a unique computational 
framework 0, 0|- As discussed above, the OFDFT is 
likely to play an important role in the emerging multi- 
scale modeling and simulation of condensed matter sys- 
tems, however the criterion to measure the validity of 
such an approach is based on the quality of the kinetic 
functional employed Q. In principle, if we had the 
exact kinetic functional the variational problem based 
on the Hohenberg-Kohn theorem would not require the 
KS orbital-based approximation anymore. In practice 
we must design kinetic functionals, which are physically 
well founded and computationally less demanding than 
the KS-based approaches. The most popular currently 
used kinetic functional consists of the Thomas-Fermi 
and Weizsacker term, which are local, plus a linear re- 
sponse term which is non local Recently one of 



the authors of the current work has explored a differ- 
ent path than that above; this is based on obtaining a 
kinetic functional for a system of N electrons whose two 
central quantities are the one particle electron density, 
p(r), and the (N — l)-conditional probability density, 
/(r 2 , rjv|ri) @, [l(J. The first is defined as: p(r) = 
p(ri) = N f nN i V>*(ri,r 2 , ...r N )il)(r 1 ,r 2 , ...r N )dr 2 ..,dr N 
(where fijv-i is the 3 (TV — 1) dimensional spatial do- 
main of the electrons and ip is the 3(iV — 1) dimensional 
wavefunction of the electronic ground state) . The sec- 
ond term represents the probability density of finding an 
N — l electron configuration, (r 2 , rjy), for a given fixed 
value of ri. The function / satisfies the following prop- 
erties (see also [Hj]): 



/(r 2 , r N \r 1 )dr 2 ...dr N = lVri (1) 



/(ri, ...r i ...r i _i,r i+ i...,r A r|rj) = 0;for i = j;Vi,j = 1,N 

(2) 

/(ri,...ri...,rj...rfc_ 1( r fc+ i...,rjv|rfc) = 0;for i = j;Vi,j ^ k 

(3) 

In terms of p and / the equation to determine the ground 
state of the systems writes [j| : 



E = min (mm(T[f,p\) + dn + J v^p^dn 

where v(ri) is the external potential and 



T[f, P ]= / P (n) 



fijV_i 



|V ri /(r 2 , ....,rjy|ri) 
/O2, ,rjv|ri) 



■dr 2 ....dr N + (N-1) 



/(r 2 , ,rjv|ri) 



Qm — i 



ri - r 2 



dr 2 ....dr 



N 



(4) 



dri. (5) 



EqU is expressed in atomic units, i.e. h, to, e are all equal to 1. This approach, up to this point is formally rigor- 
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ous; however, to proceed further it requires an explicit 
expression for /. This latter is then built on the basis of 
the mathematical prescription of Eqs lll2l3l and a physi- 
cal first principle argument. The physical argument used 
in this case is rather simple, though powerful, that is in 
contrast to the case of a gas of noninteracting electrons, 
the gas of interacting electrons is characterized by the 
presence of the Coulomb electrostatic term acting pair- 
wise between electrons. The next step consists of writing 
/ in terms of such an interaction in a way that Eqs ll|2l3l 
are satisfied [13] • The proposed analytic form, consistent 
with the requirements above, reads 0, [l(| : 

/(r 2 , ...t n \t) = n^^e 1 ^-^^ xU i>m e^ EH ^ 

(6) 

with the normalization condition: 



an optimal choice of the / by minimizing T w.r.t 7. The 
advantage w.r.t. previous methods is that in this case the 
effects of the correlations between the electrons regard- 
ing the kinetic functional are automatically incorporated 
into a local functional. Moreover, within the OFDFT ap- 
proach, the optimization procedure for / and the reduc- 
tion of the functional from non-local to local, can be also 
used as an intermediate step within a full self-consistent 
procedure as suggested by EqJU This means that the 
kinetic functional for an OFDFT algorithm can be also 
calculated on the fly without any adjustable empirical 
parameters. This, as Eq|4] shows, is valid beyond the 
case of uniform gas that here we use as an illustrative 
example 



.dr 



N 



(7) 

where Ejj^i^Tj) — ^j^z^-p and 7 is a free parameter 
to be optimized by minimizing T (w.r.t. 7) for a fixed 
p. The resulting kinetic functional consists of a local 
term (the Weizsacker term, the second term on the r.h.s. 
of EqU]) and an N— dimensional non local term known 
as the non-local Fisher Information functional, which we 
will indicate as I[p] [l^ (first term on the r.h.s. of EqJ5]). 
The other term in T is the Coulomb electron-electron 
term, which we will indicate as C[p], and is also a multi- 
dimensional integral. The aim of this work is to show, for 
a uniform gas how by employing a Monte Carlo sampling 
of the electronic configurations in space, the non local ki- 
netic term, I[p] (and in general the functional T[f, p\) can 
be reduced numerically to a local functional. Moreover, 
within the same approach, one can perform numerically 
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MC evaluation of the non local part of the Ki- 
netic Functional: At this point, the problem we must 
address to the aim of obtaining a local functional from 
Eq. [5] is that of reducing the 3iV-dimensional integrals to 
3-dimensional ones. The Metropolis MC approach has 
been shown to be ideal for this kind of problems (see 
standard textbooks in the field, e.g. 15], or Ref.[l6|]). In 
general, the choice of this kind of stochastic algorithms 
is a natural one for high-dimensional problems, since 
their efficiency increases, if compared to any other non- 
stochastic method, as the number of dimensions of the 
integral increases. In this case, the integration is done 
w.r.t. the variables r2, ...rjv for a fixed value of 17 at a 
given density p(ri). By performing the integration for a 
large number of fixed 17 values we then obtain the ex- 
pression of the functional local in 17 . In order to evaluate 
Eq. [5] within the MC approach , we need to rewrite it as 
foUows:r[/,p] - / p(rt) [J(p(ri)) + C(p(r x ))] dr u where 
the key quantities I(p(ri)) and C(p(ri)) at a fixed 17 can 
be conveniently written in the following form: 



I(p(r 1 )) + C(p(r 1 )) = 



dr N ~!f 



1 

N 



EE 

= 1,N j>i 



dr N ~ 



(8) 



P(ri),7 



r 



In the above equation the two integrand functions are 
written as a product of /, that acts here as a weighting 
function, and another function (of all the coordinates). 



This is the form suitable for the importance sampling MC 
algorithm [lij]. In fact, the integrals can be stochastically 
evaluated by using: 



I(p(v 1 )) + C(p(r 1 )) = ^ 



m=l,M 



1 


Vi/ 
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/ 



n E E' 

i=l,N j>i 



(9) 



P(ri),7 



where M is the number of randomly chosen configura- tions, or configurational space points, at which the func- 
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FIG. 1: The three panels show the values of 7 that minimize the functional J(p(n)) + C(p(ri)) at three representative densities 
(for N=100 electrons). Energies and densities are expressed in atomic units. 



tion in the square brackets is sampled. The sum in Eq. [5] 
converges to the integral of Eq. [8]when the configurations 
are generated according to /. In practice, since / is not 
known explicitly (the normalization factor is a multidi- 
mensional integral), a Markov chain is constructed with 
limiting density /. For the case of a uniform gas, the 
Metropolis MC algorithm consists of randomly displac- 
ing the present configuration of the system of N electrons 
and then accepting the new configuration, in a way that 
the density of particles in average stays constant, with 
probability (i.e. with acceptance rule): 



a (old — > new) = min 1, 



old 



(10) 



The robustness of the algorithm lies on the fact that the 
terms in the sum (Eq. O can be evaluated at each config- 
uration with relatively low computational cost. In fact, 
this is trivial for the l/|r,- — Tj\ term, while some algebra 



reduces the term 



Vi/// 



as a function of V\Eh (1*1 , r„ 



We modeled a uniform distribution of electrons by 



means of a system of N electron (with N=10, 25, 50, 
100, and 250) in a cubic box. The MC scheme works 
as follows: one electron is selected at random and a trial 
move is attempted. We adopt a trial move as a uniformly 
distributed displacement of a randomly selected electron. 
Explicitly, the acceptance rule is: 




FIG. 2: Calculated values of I(p(ri)) for N = 50, 100, and 
250 electrons. The continuous line is the fit we propose (see 
text). 



where k labels the attempted moved electron. If the move 
is accepted the terms of the sum in Eq. [9] are evaluated 
for the average in the new configuration; otherwise an- 
other instance of the quantities evaluated in the old con- 
figuration entries the average. Periodic boundary condi- 
tions and minimum image convention were imposed. In 
practice, each displaced electron was in turn the center of 
the box and only one instance of each particle was used in 
evaluating the relevant quantities in equations QT] and [9] 
Periodic replicas of the system would be necessary for the 
evaluation of the slowly decaying "Coulomb" dependence 
of C(/3(ri)); in contrast, for the evaluation of I(p(ri)), it 
would be physically wrong to count correlations of pe- 
riodic replicas (therein considering also self-correlations 



between the displaced electron and its images): each elec- 
tron should contribute once to the integral. For the pur- 
poses of this paper, we did not need to worry about the 
finite size dependence of C(p(ri)), since we have notably 
found out that both the value 7 that minimizes r[/, p] 
and /(p(ri)) do not depend on the size N of the system. 

For the case of uniform gas, once obtained the sum 
in Eq. [9l the integral r[/, p] is calculated by just mul- 
tiplying J(p(ri)) + C(p(ri)) by the number of electrons 
N. 

Numerical Results: We performed the calculation 
at several densities, selected in a way to resemble those 
of real metals. For each density we evaluated the sum 
in Eq. [5] at different values of 7, searching for the value 
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yielding the minimum in the sum. Fig. Q] shows the 
minimum found for T as a function of 7 for three dif- 
ferent densities (a) high, resembling the realistic case of 
Pt (~ 1.0 e/bohr 3 ) (b) intermediate, resembling Ti (~ 
0.2 e/bohr 3 ), (c) low, resembling Na (~ 0.04 e/bohr 3 ). 
The existence of a minimum proves that indeed the 
problem is variational and allows to find the optimal 
7 that in this way is no more an empirical parame- 
ter. Having the optimal 7 for the different densities, 
we can numerically express /(p(ri)) and then fit the 
numerical results to a functional form so that we can 
have an analytic expression. The formula we obtain is: 
J(p(ri)) = A + Bln(p(ri)), with A = 0.860 ± 0.022 and 
B = 0.224 ± 0.012. The calculated values of /(p(ri)) at 
the optimal 7 are plotted in Fig. [2] together with the fit. 
Finally, we can propose a new local kinetic functional for 
slowly varying densities, ready for OFDFT based codes: 

K[p] = \J ]Vp p { ^ 2 dv l+ J p{v 1 ){A + B\np{v 1 ))dv 1 . 

(12) 

In conclusion, we have illustrated a computational pro- 
tocol which combines analytical and stochastic numeri- 
cal approach to address the problem of designing local 
kinetic functional for OFDFT based codes. We have ap- 
plied it to the test case of uniform interacting gas of spin- 
less electrons and given arguments from which it emerges 
the general character of the approach for any system in 
the ground state. Anyway, it must be said that despite 
being an ideal system, the uniform gas has been often 
employed to develop tools for realistic calculations (see 
e.g. the very recent work of Drummond and Needs [18l]). 
often under the approximation of spinless particles. For 
this reason the specific results of this work are already 
relevant beyond the role of illustrative example. The ap- 
proach shown here is internally (mathematically) consis- 
tent and, computationally not expensive and rather flex- 
ible. It may represent a direction along which to proceed 
since it can be technically improved by both the math- 
ematical and computational point of view by improving 
the functional form of / (as discussed in Refs. 12, [Hj]) 
and the sampling of the electronic configurations. 

LMG acknowledges the AvH foundation for financial 
support. We thank M.Praprotnik for a critical reading 
of the manuscript. 



[1] Y. A. Wang and E. A. Carter, " Orbital- Free Kinetic- 
Energy Density Functional Theory," in "Theoretical 
Methods in Condensed Phase Chemistry," edited by S. D. 
Schwartz (Kluwer, Dordrecht, 2000), Chap. 5, pp. 117- 
184. 

[2] S. Watson, and E. Carter, Comput.Phys.Comm. 128, 67 
(2000). 

[3] N.Choly, E.Kaxiras, Solid State Communications 121, 
281, (2002). 



[4] J.D.Chai, J.A.Weeks, J.Phys.Chem.B 108, 6870 (2004). 

[5] V.Gavini, K.Bhattacharya, M.Ortiz, J.Mech.Phys.Sol. 
55, 697 (2007). 

[6] R.L.Hayes, G.Ho, M.Ortiz and E.A.Carter, Phil.Mag. 
86, 2343 (2006). 

[7] R.L.Hayes, M.Fago, M.Ortiz and E.A.Carter, Multiscale 
Model.Sim. 4, 359 (2005). 

[8] T.J.Frankcombe, G.J.Kroes, N.I.Choly, E.Kaxiras, 
J.Phys.Chem.B 109, 16554 (2005). 

[9] L.Delle Site, J.Phys.A:Math.Gen. 39, 3047 (2006) 
[10] L.Delle Site, J.Phys.A:Math.Theor. 40, 2787 (2007). 
[11] P.W.Ayers, J.Math.Phys. 46, 062107 (2005). 
[12] Of course, since among the properties required for / there 
is also that of mimicking the fermionic character, the in- 
direct interaction of the Pauli principle is also included. In 
case of spinless particles, the Pauli principle is included in 
a crude way, that is by the fact that two particles cannot 
have the same state or position r. In principle one may 
design an / which can separate the interaction coming 
from the Pauli principle and that from the Coulomb term. 
Here, for simplicity, we have used only the Coulomb term 
since this form of / does capture the essential physics and 
it is numerically simple. However, in principle it would 
not be difficult to include specific effects (see also [14J). 
[13] S.B. Sears, R.G.Parr, and U.Dinur, Israeli. Chem. 19, 
165 (1980). 

[14] In particular, although the hypothesis of spinless parti- 
cles is often used in many applications of DFT, it would 
not be difficult, thought computationally more demand- 
ing, to incorporate the spin explicitly in the current ap- 
proach. In this case one may follow a first principles argu- 
ment where although the Coulomb interaction between 
electrons is still the central point, however it must be ex- 
pressed taking into account the effects of the spin. This 
may be achieved by modifying / so that electrons with 
same spin tend to avoid each other (in a long range fash- 
ion) while those with opposite spin tend to "condense" 
but the electrostatic interaction keeps them apart (in a 
short range fashion) [19L I20I. I2H . Following this principle, 
/ can be written in the same way but factorizing / in a 
term considering couples with the same spin and another 
considering those with opposite spin. They may have dif- 
ferent constants, a, f3, multiplying the term in the expo- 
nential. In turn / can be optimized minimizing V w.r.t. a 
and f3. Finally, by imposing that the total spin is zero one 
can perform the sampling of the electronic configurations 
considering, for each spatial configuration, also the spin 
configurational space. 

[15] B.J.Smit and D.Frenkel, Understanding Molecular Sim- 
ulation, Academic Press Inc. U.S. (1996). 

[16] The general principles of the integration procedure are 
similar to the ones employed (for example) for variational 
quantum MC [fjj]- However there are technical aspects 
which are peculiar of our approach as discussed in the 
text. 

[17] D.Cheperly, G.V.Chester, M.H.Kalos, Phys.Rev.B 16, 
3081 (1977). 

[18] N.D. Drummond and R.J. Needs, Phys. Rev. Lett. 99, 

166401 (2007) 
[19] E.Wigner, Phys. Rev. 46, 1002 (1934). 
[20] E.Wigner and F.Seitz, Phys. Rev. 46, 509 (1934) 
[21] D.Pines, Phys.Rev. 92, 626 (1953). 



